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Regions of the ocean's thermocline unstable to salt fingering are often observed to host 
thermohaline staircases, stacks of deep well-mixed convective layers separated by thin 
stably-stratified interfaces. Decades after their discovery, however, their origin remains 
controversial. In this paper we use 3D direct numerical simulations to shed light on the 
problem. We study the evolution of an analogous double-diffusive system, starting from 
an initial statistically homogeneous fingering state and find that it spontaneously trans- 
forms into a layered state. By analysing our results in the light of the mean-field theory 
developed in Paper I, a clear picture of the sequence of events resulting in the stair- 
case formation emerges. A collective instability of homogeneous fingering convection first 
excites a field of gravity waves, with a wcll-dofincd vertical wavelength. However, the 
waves saturate early through regular but localized breaking events, and are not directly 
responsible for the formation of the staircase. Meanwhile, slower-growing, horizontally 
invariant but vertically quasi-periodic 7-modes are also excited and grow according to 
the 7— instability mcclianism. Our results suggest that the nonlinear interaction between 
these various mean-field modes of instability leads to the selection of one particular 
7— mode as the staircase progenitor. Upon reaching a critical amplitude, this progeni- 
tor overturns into a fully-formed staircase. We conclude by extending the results of our 
simulations to real oceanic parameter values, and find that the progenitor 7— mode is 
expected to grow on a timescale of a few hours, and leads to the formation of a ther- 
mohaline staircase in about one day with an initial spacing of the order of one to two 
metres. 
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1. Introduction 

Fingering convection is a small-scale mixing process driven by a well-known double- 
diffusive instability of stably stratified fiuids. When density depends on two components 
with opposite contributions to the overall stratification, a fingering instability may de- 
velop if the diffusivity of the stably stratified component is larger than the diffusivity 
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of the unstably stratified one. This commonly occurs in a large variety of natural en- 
vironments, the best studied one being the heat-salt system in the oceanic thermocline 
(Schmitt 1994 Schmitt et al. 2005). Other areas of application range from the Earth's 



atmosphere (Merceret 1977 Bois & Kubicki''2002) to magmatic melts (Tait & Jaupart 



1989J and the i nteriors of giant plan ets (|Guillot|,1999| ) and stars ( jVauclair||2004| |Char- 
bonnel fc Zalin|[2007| [sFancliffe et aL|[2007^ 



Scientific interest in the dynamics of fingering convection, the saturated state of the 
fingering instability, is twofold. When viewed as a small-scale turbulent process, it is an 
intrinsic source of diapycnal mixing in stably stratified regions, which operates even in 
the absence of mechanical forcing. As such, it can play an important role in controlling 
the global transport of heat and chemical species in the system considered (see reviews 
by Schmitt, 1995 and Ruddick & Gargett (2003)). Fingering convection is also known 
to drive dynamics on much larger scales by exciting a variety of secondary instabilities, 



such as internal gravity waves through the collective instability (jStern,)1969 Stern et al. 



2001[ ), thermohaline intrusions ( |Stern|1967[ [Toole fc Georgi|1981[|Walsh fc Ruddick|1995[ ) 

and the more recently discovered 7-instability fRadko 2003). Generally speaking, these 
secondary instabilities are driven by a positive feedback mechanism between large-scale 
temperature and salinity perturbations, and the turbulent fingering fluxes induced by 
these perturbations. In Paper I, we provided definitive measurements of the small-scale 
turbulent fluxes induced by fingering convection in the heat-salt system, and unified 
all previous theories of the large-scale secondary instabilities in a single "mean-field" 
framework. In this paper, we now turn to one of the longest standing unsolved problems 
related to fingering convection, namely the formation of thermohaline staircases. 

Thermohaline staircases are formed as stacks of deep convective layers, both thermally 
and conipositionally well-mixed, separated by thin fingering interfaces. They are ubiq- 
uitously associated with active fingering convection, and have been observed both in 
laboratory experiments ( Stern fc Turner||1969 Krishnamurti|[2003 2009) and in oceanic 
field measurements ( Tait fc Howe|197l Schmitt et aL|2005 ). Oceanic staircases are long- 
lived, and can span the entire thermocline, with individual layer heights ranging from 
lOm-lOOm and horizontal extents of km-size or more. Regions of the upper ocean which 
exhibit a staircase stratification have much larger diapycnal mixing rates compared with 



smoother regions supporting the same overall temperature and salinity contrasts ( Schmitt 



et al. 2005 Veronis 2007). The layering phenomenon therefore needs to be understood 
to quantify its role in the overall transport balance in the ocean, and by extension, in 
many other systems as well. 

Forty years after their original discovery, however, the mechanism responsible for the 
formation of thermohaline staircases essentially remains an enigma. Existing theories 
might be roughly classified into two groups. In a first class of models, layer formation is 
triggered by external effects in addition to fingering convection. It has been speculated 
for example that staircases and homogenous fingering both represent distinct metastable 
equilibria of the same system, and that layers are formed by external events which cause 
vigorous mechanical mixing and drive the system from one equilibrium to the other 
(Stern & Turner 19691. Since staircases are often observed in regions with significant 
lateral gradients of heat and salinity (such as the Mediterranean outflow for example), 
it has also been proposed that these gradients are actually required to trigger their 
formation, through the excitation and nonlinear evolution of the thermohaline intrusions 
dWalsh fc Ruddick|[l995l|Merryfield|[2000l ). 

By contrast, a second class of ideas favour the notion that layering is inherent to 
the dynamics of fingering convection, with no need for external forcing. In these mod- 
els, staircase formation is thought to arise from the nonlinear development of the two 
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other large-scale modes of instability mentioned earlier, namely the collective instability 
( |Stern||1969| [Stern a/.] [2001^ or the 7-instability (Radko,,2003| . [Stern et al] ( [2001[ ) 
studied the linear and nonlinear evolution of gravity waves excited by the collective in- 
stability, arguing that the overstable waves would eventually overturn and break into 
layers. However, while their numerical simulations show evidence for localized breaking 
events, these do not result in layer formation. Radko (2003) on the other hand put for- 
ward the 7-instability as the direct precursor of thermohaline staircases, since by nature 
the normal modes of this instability are vertically modulated but horizontally invariant 
perturbations in temperature and salinity. His two-dimensional direct numerical simu- 
lations (DNS), which are the first to exhibit unseeded spontaneous staircase formation 
from an initially homogeneous fingering field, appear to confirm his hypothesis. 

Despite recent progress in this type of models, a number of fundamental questions 
remain to be addressed. The first relates to the validity of the mean-field framework: 
do the mean-field equations correctly model the large-scale dynamics of fingering sys- 
tems? This question has so far only partially been answered for the linearised mean-field 
equations through idealised tests in which a single large-scale perturbation is seeded on 
the fingering field, its growth rate monitored and successfully compared with theoretical 
predictions (Stern et al. 2001 Radko 20031. However, whether such perturbations are 
indeed excited by the turbulence with a sufficiently large amplitude to rise above the 
noise level, and with sufficient scale separation for mean-field theory to be applicable, 
remains to be determined in the general case. Furthermore, we note that these tests 
have so far only been performed in two-dimensional numerical simulations. Given the 
well-known pathological nature of the energy transport from small to large scales in 2D 
turbulence ( Kraichnan 1967 ) , convincing validation of mean-field theory can only come 
from 3D simulations. 

The second question is related to the problem of layer formation itself. Existing theo- 
ries ( Stern et aL|2001 Radko[[2003 1 have each put forward a particular mean-field mode 
of instability as the precursor for staircase formation. As such, they consider each mech- 
anism in isolation and ignore other possible unstable modes. In Paper I we identified the 
fastest growing modes, and showed that in the parameter regime where staircases are 
observed several kinds of mean field instabilities can simultaneously be excited. In fact 
modes on a considerable range of spatial scales, both direct and oscillatory in nature, are 
predicted to grow by mean field theory, and it seems likely that their mutual interaction 
disrupts their evolution. The resulting dynamics remains to be explored. 

Finally, all previously proposed models assume that the nonlinear evolution of a mean- 
field mode eventually triggers the formation of layers. Unfortunately, the predictive power 
of these mean-field models breaks down as soon as turbulence on a wide and continuous 
range of spatial scales is generated. Three-dimensional computer simulations resolving 
all relevant scales currently are therefore the only way to study the final stages of layer 
formation. 

The purpose of the present paper is to answer some of these questions, and ultimately 
illuminate the problem of staircase formation, via large-scale, three-dimensional numer- 
ical simulations of fingering convection. Since both small-scale fingers as well as deep 
layers have to be resolved simultaneously for this to be possible, and since fingers and 
staircases evolve on vastly disparate timescales, such simulations had until recently been 
impossible. Thanks to the rapid advances in high-performance computing, this is no 
longer the case. 

Our paper is organised as follows. We begin by briefly summarising in Section [2| the 
basic equations for fingering convection, and describe our numerical model setup and 
selection of governing parameters. In Section[3|we present our results, the first large-scale 
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three-dimensional simulation showing spontaneous layer formation, which is fully resolved 
down to the dissipative scale. We analyse the simulation in the light of the mean-field 
theory presented in Paper I, focusing on the evolution of the various large-scale modes 
of linear instability. Our results enable us to place constraints on the applicability of 
mean-field theory, and lead to a clear and simple picture of the sequence of events finally 
resulting in layer formation. We conclude in Section |4] by summarising our results and 
by discussing their implications both in the oceanographic context, and further afield. 



2. Equations and model setup 

2.1. Governing equations 

In natural systems, both fingering convection and staircase formation usually occur far 
from physical boundaries. It is therefore important to select a numerical setup which 



minimises boundary-effects. Following Stern et al. (2001) and Radko (2003), for exam 



pie, we consider a fingering-unstable system, permanently forced by background vertical 
temperature and salinity gradients Tqz and Sqz- In this case, the non-dimensional model 
equations are: 

^ (^l^+u-Vu^ =-Vp+(T-5)k + V2u , (2.1a) 
V-u = 0, (2.16) 

dT 

— + W + U-VT = V^T , (2.1c) 

f)Q 1 

-w + u-WS =^tW^S , (2.1d) 



dt R, 







where T denotes perturbations away from the linearly stratified background temper- 
ature, S the similarly defined salinity perturbation, p the pressure perturbation from 
hydrostatic balance, and the vector u = (u,v,w) is the velocity vector. The above sys- 
tem has been non-dimensionalised using the anticipated finger scale d = {ktv / gaTozY/^ 



(Stern 1960) as the lengthscale (where kt is the thermal diffusivity, u the kinematic 
viscosity, g is gravity and a denotes the thermal expansion coefficient), and the thermal 
diffusion timescale across d, namely d^ /kx, as the timescale. Temperature and salinity 
have been rescaled by T^zd and aT^zd/ {3 respectively, where j3 is the saline contraction 
coefficient. Three non-dimensional parameters control the behaviour of the system: the 
diffusivity ratio r = ks/kt, where ks is the salt diffusivity, the Prandtl number Pr 
= v/kt, and the background density ratio i?o — ciToz/ PSoz- 

To minimise the influence of boundaries, perturbations are assumed to be periodic in 
all three spatial dimensions. Hence for a computational domain of size (L^, Ly, Lz) we 
set q{x, y, z, t) = q{x + L^;, y, z, t) = y + Ly,z,t) = q(x, y,z + Lz,t) for q G {T, 5, u}. 
This configuration guarantees in parti cular th at layers cannot be generated artificially by 



flux convergence at soHd boundaries ( Ozgokmen et a/.||l998[ |Merryfield||2000| . It is also 



important to note that the total horizontal averages of the temperature and salinity fields 
are not assumed nor forced to be zero, so that the mean vertical profiles of temperature 
and salinity can evolve freely with time. 

Finally, note that by contrast with Paper I, we do not include lateral gradients, and 
thus suppress the mechanism driving intrusive modes. As such we restrict our analysis to 
spontaneously emerging staircases, and show how they may form even in the absence of 
any additional forcing. Although intrusions are important in regions of the ocean subject 
to very strong lateral gradients, this assumption is reasonable in more typical regions of 
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Figure 1. Numerically determined buoyancy fluxes due to heat and salt, Ft = (wT) and 
Fs = (wS) respectively, their ratio 7, and the Stern number A — {Fs — Ft) / Pr (1 — l/Rp), as 
a function of Ro for t — 1/3 and Pr = 7. The theoretical prediction of |Schmitt| (|1979|) is shown 
in b ) for reference. The fluxes were measured after integrating the model described in Section 
2.1|to saturatio n, in a small domain of size 5 x 5 x 15 fastest-growing finger widths (FGWs, see 
Schmitt (1979 1 ). The experimental protocol for these measurements is simil ar to that d escribed 
in the Appendix of Paper I. A necessary condition for the 7-instability _|Radko||200"3l ) is that 
7 = Ft /Fs should be a decreasing function of Rq. As shown in Figure 1 a, this ha ppens here 
when 7?o < 1-7. A necessary condition for the collective instability ( [Stern "ei aZ.|200l[ ) is that the 
Stern number A should be larger than an order one factor. As seen in Figure Ofc this condition 
is satisfied only for Rq < 1.7 as well. 



the thermocline since Paper I shows that for smaU-to-moderate lateral gradients, intrusive 
modes grow on much slower timescales than either gravity-wave modes or 7— modes. 

2.2. Selection of the governing parameters 

Simulations which have to resolve both fingers and possible larger scale structures such 
as thermohaline staircases still pose significant numerical challenges. As seen in Paper 
I (see Table 1), for parameter values appropriate for salty water (Pr = 7, t ~ 0.01), 
even a "small domain simulation" containing only a few fingers requires a resolution of 
the order of 1000"^ for low values of the density ratio to be fully resolved. Since our goal 
is to simulate a much larger domain, we are constrained to use a larger value of the 
diffusivity ratio instead, and choose t = 1/3. The Prandtl number on the other hand can 
be taken to be that of water (Pr = 7) without difficulty. Guided by Radko (2003), we 
adopt i?o — 1.1 as the background density ratio. 

As described in the introduction, we are interested in using the numerical simulation 
to study the causal relationship between the collective- and 7-instabilities and layer 
formation. Oceanic staircases are observed in regions with low density ratio, typically 
Rq < 1.8 (e.g. Schmitt 1981). In this parameter regime, mean-field theory predicts that 
gravity- wave modes grow much more rapidly than 7-modes (see Figure 5 Paper I). We 
must then verify that our selected numerical parameter values Pr = 7, r = 1/3 and 
i?o — 1-1 yield the same relative ordering of the growth rates of the various mean- field 
modes of instability. As shown in Paper I, these growth rates are solutions of a cubic 
equation where the coefficients of the cubic depend on the governing parameters {Rq, 
Pr, r), on the spatial structure of the mode, and on the functional dependence of the 
turbulent heat and salt fluxes on the density ratio. Using the methodology outlined in 
Paper I, we first measure these flux laws from small-domain experiments (Figure [T]), and 
then use the results to calculate the growth rates of the mean-field modes (see Figure [2]). 

Figure [2] confirms that gravity-wave modes and 7-modes are both unstable in our 
numerical parameter regime, providing us with the opportunity to study them simul- 
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Figure 2. Growth rate A of mean-field instabilities for Pr = 7, r = 1/3, and _Ro ~ 1-1, using the 
mean-field theory presented in Paper I and the flux laws presented in Figure [l] The horizontal 
axis represents the vertical wavelength in units of d, and the different curves represent different 
horizontal wavelengths as noted in the caption. Note that mean-field theory does not apply to 
structures which only contain a few individual fingers. This limits its val idity to modes with 
horizontal and/or vertical wavelengths larger than about lOOd (see Section 3.1 for detail). 



taneously. Furthermore gravity waves are expected to dominate the system's behaviour 
while the 7-modes grow much more slowly. This choice of parameters thus guarantees 
that the basic properties of the unstable modes are comparable to those of the heat-salt 
system. Hence, despite the difference in parameters, we anticipate that our simulation 
should provide a reasonable view of the dynamics of real thermohaline staircase formation 
at Pr = 7, T = 1/100 and Rq - 1.5. 

2.3. Numerical considerations 

Our large-scale simulation is performed in a domain of size 335d x 335d x 536d, which cor- 
responds to 27 X 27 X 43FGWs. The numerical algorithm for the solution of ( 2.1a|2.1d ) 



with periodic boundary conditions is based on the standard Patterson-Orzag spectral 
algorithm and is presented in more detail in Paper I. Using small-domain tests to de- 
termine an adequate resolution for this parameter regime, we settled on a resolution of 
576'^ equivalent grid points for the large-domain run. The computation was initialised 
with small amplitude white noise perturbations in the temperature and salinity fields, 
and evolved from t — to t — 900 (in dimensionless time units). Note that the initial 
resolution was sufficient for most of the simulation but became inadequate in the layered 
regime, leading to a considerable accumulation of energy in the high vertical wave-number 
modes of the salinity field. Doubling the vertical resolution before layers start to form 
solved this problem: the results were interpolated to a 576 x 576 x 1152 grid slightly 
before the layering transition and the simulation was then re-computed on the finer grid 
from there on. 



3. Direct simulation of staircase formation 

The dynamics of the system, as observed in our numerical simulation, are divided 
into three distinct phases illustrated in Figure [3j Note that a movie of this simulation 
is available in the online supplementary material. In Phase I, initial perturbations are 
amplified by the fingering instability, grow and eventually saturate into a state of vigorous 
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Figure 3. (a-c) Simulated layer formation in a large triply-periodic domain (335dx 335dx 536d). 
Temperature perturbations normalised to the temperature difference across the whole domain 
are shown in colour, with their respective colour scales in each panel. The top panels are visual- 
isations on the data-cube faces, while the lower ones are volume-rendered images. Snapshots are 
shown for three characteristic dynamical phases. Fingering convection develops first (Phase I) 
and soon becomes unstable to large-scale gravity waves saturating at finite amplitude (Phase II). 
Eventually, vigorously convecting layers form, separated by thin fingering interfaces (Phase III), 
(d) Time-series of the Nusselt number Nu= 1 — Ft and of the turbulent fiux ratio 7 = Ft/Fs- 



fingering convection (0 ^ t ^ 100; Figure [3^). The global heat and salt fluxes, measured 
by Nu — 1~Ft and 7 = Ft/Fs, are found to be identical to those determined previously 
in a smaller computational domain (see Figure[l^). However, the system does not remain 
in this homogeneous fingering state for long. Gravity waves rapidly emerge, grow and 
saturate (Phase II, for 100 ^ t ^ 600), later followed by a sharp transition to a layered 
state (Phase III, t > 600). Phases II and III are now analysed in more detail. 

3.1. Phase II : gravity-wave phase 

Shortly after saturation of the fingering instability (t ^ 100), large-scale oscillatory modu- 
lations of the fingering field appear. They rapidly grow until t = 150 when they, in turn, 
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also saturate. We argue below that these modulations are large-scale finite amplitude 
gravity waves excited by the collective instability. 

The visualisation of the temperature field presented in Figure [3]d and in the online 
movie, for t > 150, reveals the presence of a dominant perturbation pattern, with a hor- 
izontal wavelength commensurate with the box width (335(i), and a vertical wavelength 
half the box height (268(i). Its temperature amplitude at saturation is of the order of 
15% of the background temperature difference across the domain, and equivalently 30% 
of the temperature difference across one wavelength. The spatial pattern and relative 
amplitudes of the density and salinity perturbations are similar. 

Figure [3]i reveals that these large-scale perturbations strongly affect the global heat 
and salt fluxes. Both Nu and 7 begin to show visible oscillations around their respective 
background values (i.e. as achieved in the homogeneous fingering state) shortly after 
t = 100, with increasingly larger amplitudes until saturation at about t — 150. By analogy 
with the effect of a pure gravity wave on global heat transport, we expect the Nusselt 
number to oscillate at twice the wave frequency, and similarly for 7. A power spectrum 
of the Nusselt number reveals a strong peak at w ~ 1.04, so we estimate the wave 
frequency to be ujg ~ 0.52. In our dimensionless units, the theoretical oscillation frequency 
for a pure gravity wave with phase planes inclined at an angle 6 from the vertical is 
Ljg = •^Pr(l — I/Rq) cosO, in other words Wg = 0.8 cos for our selected parameter 
values. The dominant frequency observed is therefore compatible with the presence of a 
gravity wave with 9 — 49°, which is consistent with our qualitative observation of the 
wave geometry. These results also confirm the theoretical prediction that gravity waves 
are the most rapidly growing mean-field mode in this parameter regime, and suggest 
that their absence from the 2D simulations of Radko (2003) is an artefact of reduced 
dimensionalitjif] 

We now study the gravity waves more quantitatively in the light of the linear mean-field 
stability analysis presented in Paper I. In this triply-periodic system, all mean-field modes 
have the normal form qimk exp(ilx + imy +ikz) where qimk is the mode amplitude, so that 
their predicted evolution can be compared straightforwardly with Fourier-transforms of 
the output of the simulation. Using this method, we extract from the numerical results 
the amplitude of each velocity component u, v and w, and compute the kinetic energy 
Eimk = (l^/mfeP + IwimfcP + |w/mfep)/2 of any mode with wave-vector k = {I, m, fc), as a 
function of time. The results are shown in Figure [4] 

For simplicity, in what follows, we use the notation Z„ = 2mT/Lx, rrin — 2m: / Ly, as well 
as kn = 2m:/ Lz. Hence the dominant gravity- wave modes identified by eye in Figure [sjj, 
which have one wavelength in the horizontal direction and two in the vertical direction, 
are denoted by k = (Zi,0, ^2) or k = (0,TOi,fc2). Figure [4] shows that this set of modes, 
identified by a star symbol, begins to grow exponentially as soon as the background 
fingering convection reaches a saturated state. By t = 200, their kinetic energy is ten 
times larger than any of the other modes. Comparison of the observed growth rate with 
mean-field theory shows reasonable agreement from t = 50 to t = 150, although is slower 
than predicted by about 30%. 

It is interesting to note that this dominant set of modes is not the most rapidly growing 
one according to linear theory: as seen in Figure |4] a variety of other gravity- wave modes 
have equal or larger predicted growth rates, as for example k — {l2,0,k^). However, 

f The difference in the wave dynamics in 2D and 3D may be attributed to the stronger 
negative feedback of the wave shear on salt fingers in 2D. In three dimensions, salt fingers tend 
to align into salt sheets parallel to the shear flow, which only moderately affects their vertical 
transport. In 2D such alignment is geometrically impossible and therefore adverse effects of 
shear on salt-fingers, and thus on the collective instability, are more profound. 
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Figure 4. Kinetic energy in the first nine quartets of gravity wave modes, compared witfi tlie 
growth rate predicted using the mean-field theory presented in Paper I. Owing to the horizontal 
symmetry of the problem, modes (±i, ±m, k) and (±m, ±1, k) all have the same predicted growth 
rate, and are superimposed in each panel. The vertical wavenumber increases from top to bottom 
and the horizontal wavenumber increases from left to right. The (±Zi,0,fc2) and (0,±mi,A;2) 
modes (denoted by a star) are the "dominant modes" referred to in the main text. 



not all modes grow as predicted. A systematic comparison between the observed and 
predicted growth rates of all nine sets of modes shown in Figure |4] reveals that mean- field 
theory becomes progressively less accurate with increasing wavenumber. This information 
is summarised in Figure [5] Since for our selected set of parameters a finger is typically 1 
FGW wide and 5 FGW high (see Paper I), and IFGW ~ 13d, our results suggest that 
modes grow as predicted only if their horizontal extent is wider than about 10 fingers and 
their vertical extent is larger than about 5 fingers. Modes with intermediate horizontal or 
vertical extent still grow but slower than expected. Modes which only contain 3 or fewer 
fingers (horizontally or vertically) do not grow at all. These results are not surprising 
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Figure 5. Contour plot of the growth rate of mean- field modes for a) our selected parameters 
Pr = 7, r = 1/3, i?o = 1.1 and b) the heat/salt case Pr = 7,r = 1/100 with Ro = 1.5. The 
individual modes presented in Figure|4]are identified on the left plot by symbols, showing modes 
which grow with the expected growth rates as filled dots, modes which grow but more slowly 
than expected as empty dots, and modes which do not grow at all as crosses. The individual 
modes, from top to bottom, have k = ki, k = k2, k = k^, and k = ki, and from right to left, 
have {I = h,m = 0), [l = h,m = mi), {I = Z2, m = 0), [I = h,m = 0) and (Z = Z4, m = 0). The 
horizontal and vertical lines show the size of the structures considered in terms of finger widths 
and heights (taking, as shown in Paper I, that fingers are about 1 FGW wide and 5 FGW tall). 
The equivalent plot for the heat/salt system shown on the right, which was obtained by using 
the turbulent fiux laws presented in Paper I, is very similar when viewed in terms of finger sizes. 



given that separation between the mode scale and the individual finger scale is required 
for mean-field theory to be applicable. Note that, when viewed in terms of finger sizes, the 
growth rate plot for the heat /salt case at Rq = 1.5, shown in figure [sjD), is very similar to 
the one obtained for r = 1/3. It therefore seems reasonable to expect that finger driven 
oceanic gravity waves will also occur on length scales larger than the predicted fastest 
growing mean-field modes. 

Finally, Figure |4] also reveals that all gravity- wave modes stop growing roughly at the 
same time, around t ~ 150 — 200. We now show that this is due to localised overturning 
events which are regularly triggered when some of the dominant modes constructively 
interfere to produce very high amplitude perturbations in the density field. 

A single gravity-wave mode with vertical wavenumber k causes quasi-periodic oscil- 
lations of the total density field (background -I- perturbations). As the mode amplitude 
Pi.m.k grows, localised regions are progressively more weakly stratified. They become un- 
stable to direct overturning convection when the perturbation reaches the critical value 
Pi,m,k — (1 — l/-Ro)/^. Figure[6^ compares the spectral power in the density perturbation 
of one of the four dominants modes - defined as the norm of the relevant coefficient of 
the Fourier transform of the density field - with the square of the critical amplitude 
for overturning. For this single mode, the criterion for overturn is clearly never reached. 
However, four of these dominant modes co-exist and can interact constructively to cause 
a localised overturning event. This effect is demonstrated in Figure [6]3, which shows the 
filling factor of large-scale regions with positive density gradient, or in other words the 
relative volume of regions unstable to overturning convection. As expected, this filling 
factor is close to zero in the homogeneous fingering phase. It then rapidly increases to 
a few percent around t — 150, which corresponds to the time when the density spectral 
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Figure 6. Top: Temporal evolution of one of the dominant gravity-wave modes ((Zi, 0, ^2), with 
horizontal wavelength 335d, vertical wavelength 236d) and of the (0,0,^2) 7-mode (vertical 
wavelength: 236d) respectively, as simulated numerically and as predicted by mean-field theory 
(Paper I). Shown is the norm of the Fourier coefficient of the density perturbation (the density 
spectral power). While the gravity- wave mode rapidly saturates the 7-mode slowly but steadily 
grows up to the critical amplitude for overturning (black line). Bottom: Filling factor of convec- 
tively unstable regions {dp/dz > 0), and total kinetic energy of the system. Note that in order 
to eliminate finger-scale density inversions, the filling factor is obtained from a low-pass filtered 
density field. 

power of one of the dominant gravity-wave modes reaches about a tenth of the critical 
amplitude for overturning. After t = 150, the steady conversion of kinetic energy into 
potential energy by the localised and intermittent overturning events causes all the 
gravity- wave modes to stop growing and not just the dominant ones (see Figure |4|. This 
result is not surprising given the nonlinear nature of the breaking events. 

Our results show that the gravity wave-field saturates at fairly high amplitude, but 
without directly triggering the formation of thermohaline layers. The gravity- wave-dominated 
phase continues on for many more oscillation periods (from t ~ 200 to t ^ 600) and the 
movie of the simulations (see online material) shows no obvious visible change in the 
statistical properties of the flows. Our simulations show that the collective instability, 
as was originally proposed, does not appear to be responsible for the transition to the 
layered state observed in Phase III. 

3.2. Phase III : Layered phase 

The gravity-waves dominated phase suddenly ends around t ~ 600, with the spontaneous 
emergence of two convecting layers separated by thin fingering interfaces (Figure [s];) . 
This layering transition is accompanied by a significant increase in the global heat and 
salinity fluxes (Figure [Sji) , as well as in the total kinetic energy (Figure [6]). Figure [T] 
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Figure 7. Horizontally averaged total density and salinity (in both cases, background + per- 
turbation) and standard deviation of salinity cr{S). Density and salinity are normalised by the 
respective background contrast over the entire domain, and height is given in units of the domain 
height. The three lines correspond to the visualisations in Figure [3] In the fingering phase, the 
salinity and density profiles are indistinguishable from the background profiles, and the standard 
deviation of salinity is statistically homogeneous. In the gravity-wave phase, the profiles remain 
close to the background profiles, but the presence of a (0, 0, ^4) 7— mode can be seen in the o-(S') 
profile. The layered phase exhibits a clear staircase in the salinity and density profiles, while 
ct{S) shows strong peaks corresponding to the position of the fingering interfaces. 

shows the resulting staircase in the mean profiles: turbulent mixing leads to a nearly 
homogeneous temperature and composition in the bulk of the layers, while the strongly- 
stratified interfaces exhibit larger spatio-temporal fluctuations in the temperature and 
salinity (as measured by the standard deviation from the mean) caused by the remaining 
fingering. The staircase, once formed, is extremely robust and persists despite being 
occasionally pierced by strong convective plumes. 

3.3. The 'y -instability and the formation of thermohaline staircases 

In order to understand what causes the layering transition, we now study Phase II again 
in the fight of predictions from mean-field theory. This time, we focus on the 7-modes 
discussed in Paper I and in Section [T] By virtue of being horizontally invariant and 
vertically periodic, 7-modes might straightforwardly overturn into a fully-formed ther- 
mohaline staircase once they reach the same critical amplitude as the one discussed in 
the context of gravity- wave modes above (see Section [O] ). 

Since the emergent staircase in our 3D simulation has two layers, we naturally begin by 
studying the (0, 0, fc2) 7-mode. Figure |6] shows that the (0, 0, mode is indeed present 
in the simulation, and began to grow at the same time (f ~ 50 — 100) as the gravity 
waves studied above. Remarkably, we find that its growth rate is well modelled by linear 
mean-field theory all the way to about t — 600, even though the mean field approach 
does not account for the turbulent mixing induced by the overturning waves. This result 
establishes the predictive power of our formalism. 

Figure |6] shows that this 7— mode reaches the critical amplitude for the onset of convec- 
tive overturning at t ~ 580 and stops growing shortly thereafter when a regular staircase 
in the density (or temperature and salinity) profile appears (t 600) . The layers rapidly 
become fully convective, as can be seen in the strong increase in the filling factor of 
convectively unstable regions and of the total kinetic energy (see Figure [6]), and the cor- 
responding increase of the Nusselt number (see Figure [Sjj). Finally, note that the gravity 
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Figure 8. Density spectral power in the four largest vertical modes (with respective vertical 
wavenumbers ^1,^2,^3 and k4), compared with mean-field theory. The two gravest modes (top 
panel) continue growing throughout the gravity-wave-dominated phase, while the other two 
modes (lower panel) stop growing around t ~ 150 — 200. 



waves can no longer exist in the layered state. Their amplitudes are strongly reduced, 
and stop showing a clear oscillatory signal. 

3.4. The interaction of j— modes and gravity-wave modes 

The formation of thermohaline staircases in our simulation thus appears to depend only 
on the growth and eventual convective overturning of 7— modes, in a way which can be 
understood quantitatively using the mean-field theory and small-scale flux laws presented 
in Paper I and in Figure 1. In many ways, the fact that mean-fleld theory holds for the 
staircase progenitor, the (0, 0, fc2) 7— mode, is rather remarkable given the presence of the 
strong gravity wave field, and deserves further attention. Moreover, as shown by Radko 
(2003), the 7-instability suffers from an ultra-violet catastrophe whereby the mode growth 
rates increase quadratically with wavenumber. As a result, the (0, 0, fc2) 7— mode is not 
the most rapidly growing one according to linear theory, and one may wonder what led 
to its selection as the staircase progenitor. 

To answer these questions, we begin by extracting information on other 7— modes 
in the system, with vertical wavenumbers ranging from ki to (modes with higher 
wavenumbers, as discussed earlier, are poorly represented by mean-field theory). Their 
temporal evolution is shown in Figure [8] and compared with theoretical expectations. 
As expected, the (0,0, fc4) mode grows the fastest among the ones studied, followed by 
(0, 0, ^3), (0, 0, fc2) and finally (0, 0, fci). At early times [t < 200), we see that the measured 
growth rates are indeed close to the corresponding theoretical ones for all modes except 
for (0, 0, fci), which unexpectedly grows significantly faster than predicted - the origin of 
this discrepancy remains a puzzle. 
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Figure 9. Horizontally averaged density perturbation {p)h{z,t) as a function of time from 
t = to t = 800. Four intervals are shown, each on a different colour scale to provide adequate 
contrast. From t = Otot= 100, mean perturbations to the background density are close to zero. 
However, notably regular oscillations appear around t = 100, with a k4 vertical wavenumber and 
twice the temporal frequency of the dominant gravity-wave mode. By t = 200, the system shows 
clear evidence for the presence of a (0,0,^4) 7— mode, strongly modulated by the gravity-wave 
field. Around t — 500, the (0, 0, ^2) mode emerges as the dominant mode, continues to grow and 
begins to overturn around t — 600 into two mixed layers separated by sharp interfaces. 



The subsequent evolution of the modes, beyond t = 200 in FigurejS] reveals an interest- 
ing dichotomy between the two larger-scale modes (top panel) which both grow with the 
expected growth rate, and the two smaller-scale modes (lower panel) which stop growing 
at about t — 150 — 200, when the gravity waves saturate. The amplitude of the (0, 0, fc4) 
mode remains the largest one until about t = 500, even though it stopped growing at 
early times. Had it continued growing only slightly further, this mode would have been 
the one to cause the staircase formation, resulting in an initial layer height half of the 
one observed in the simulation. The evolution of the horizontally averaged density per- 
turbation (denoted as {p)h{z,t) from here on) shown in Figure [9] confirms more visually 
the dominance of the (0, 0, fc4) mode until about t — 500, and the transition to a (0, 0, 
mode thereafter (with a hint of (0,0, /ci) inducing a small asymmetry between the two 
layers). 

The results of Figure [8] naturally prompt us to wonder why the smaller-scale modes 
stop growing during the gravity-wave phase while the larger-scale, slower growing modes 
do not. This important question impinges on the problem of the selection of the 7— mode 
which eventually becomes the staircase progenitor, and thus determines the initial layer 
spacing. In order to proceed further and understand what may cause the smaller-scale 
modes to stop growing, we need to gain further insight on the nonlinear dynamics of the 
various mean-field modes. 

We can identify the dominant nonlinear mechanism for interaction between the various 
mean-field modes with a closer inspection of the {p)h{z, t) profile. Figure [9] reveals notably 
regular oscillations from t — 100 to t — 600 which can be traced back to the nonlinear 
buoyancy transport induced by the dominant gravity-wave modes studied in Section 
3.1 Indeed, a gravity- wave with vertical wavenumber k and frequency uig induces a 
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horizontally averaged vertical heat (or salt) flux (wT) or (wS) proportional to eKp{2ikz + 
2iujgt). Since the dominant gravity- wave mode has a k2 vertical wavenumber, the induced 
perturbations in {p)h are expected to have a fc4 wavenumber and to oscillate at twice 
the wave frequency. This oscillation is indeed observed, and is most regular before the 
wave-breaking events begin, between t = 100 and t = 150 in Figure [9] 

By virtue of being a quadratic quantity in the wave amplitude, the wave-induced 
flux is neglected in any standard linear mean-fleld theory. However, Figure |9] shows 
that it remains important throughout the entire wave-dominated phase from t = 100 
to t = 600. In particular, it is responsible for the oscillations which are seen to modulate 
the (0, 0, ki) 7— mode rather strongly, either enhancing it or suppressing it depending 
on the phase of the wave. One may therefore conjecture that the (0, 0, k^) (and presum- 
ably the (0, 0, fcs)) 7— modes stop growing as a result of nonlinear interactions with the 
dominant gravity-wave modes through the wave-induced flux. Unfortunately, a complete 
quantitative understanding of the numerical results from our simulations, in particular 
in view of answering the question raised above, can only be gained from a thorough non- 
linear analysis of our mean- field equations, a problem which is beyond the scope of the 
current paper. However, given that the dominant gravity-wave mode induces horizontally 
averaged fluxes with exactly the same vertical structure as that of the (0, 0, fc4) 7— mode, 
the conjecture seems reasonable. 

The question then remains of why the (0,0,^2) and (0,0, fci) 7— modes on the other 
hand still grow at the predicted rate despite the effect of the wave-induced flux. This 
result is even more remarkable for the (0,0, A;2) mode, which has an amplitude several 
orders of magnitude smaller than that of the dominant gravity-wave modes at early 
times. Again, without a full analysis of the nonlinear mean-field equations we are left 
to speculate on the matter. A plausible answer lies in the spatial mismatch between the 
small-scale flux modulation induced by the waves (working on the k^^ lengthscale) and 
the intrinsically larger scale of the two gravest 7— modes. However, showing this would 
again require a weakly nonlinear analysis of the mean-field equations, and is deferred to 
future work. 

The results of the simulations, however, stand beyond the speculations. Whatever 
the cause for this observed effect, it appears that the dominant gravity waves act as a 
"filter" for any 7-mode with significantly smaller vertical wavelength. If we generalize 
this result to other parameter regimes, the progenitor of the staircase is always likely to 
be a 7— mode with wavenumber commensurate with that of the dominant gravity waves 
excited by the collective instability. As such, it is interesting to note that both sets of 
modes play a role in the eventual formation of a thermohaline staircase, although in a 
way which was neither anticipated by Stern (19691, nor by Radko (20031. 



4. Discussion and prospects 

Our simulation, when combined with the results of mean-field theory derived in Paper 
I, presents a clear and simple view of the sequence of events which may lead to the 
spontaneous formation of oceanic staircases (in the absence of lateral gradients or other 
sources of mechanical mixing), and which can also be straightforwardly generalised to 
other double-diffusive systems. 

Let us first consider a homogeneous heat-salt system, with constant background tem- 
perature and salinity gradients unstable to fingering convection. At low enough density 
ratio (i?o < 4, see Paper I), the turbulent flux ratio 7 induced by fingering convection 
decreases most rapidly with Rq. This regime is simultaneously unstable to two types 
of secondary large-scale "mean-field" instabilities: the collective instability which ex- 
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cites gravity waves ('Stem"'1969 Stem et aL||2001 ), and the 7— instability which leads to 
the growth of horizontally-invariant perturbations in the temperature and salinity field 
( |Radko|[2003| . 

As shown in Paper I, mean-field theory predicts that gravity-wave modes grow more 
rapidly than 7-modes for values of Rq characteristic of regions of the ocean where ther- 
mohaline staircases are observed (i.e. 1 < Rq < 1.8). By comparison with numerical 
simulations at Pr = 7 and r = 1/3 we also found that the theoretically fastest-growing 
gravity-wave modes are actually too small for mean-field theory to be applicable. In- 
stead, the emerging dominant gravity-wave modes are slightly more extended and have 
smaller theoretical growth rates. If we extrapolate the results of Figure [5^ to parameter 
regimes relevant of the heat-salt system (seejsjD), taking Rq — 1.5 as a typical value for 
staircase regions, we find that the overall picture is very similar: modes which are large 
enough to grow with the theoretically predicted rate are very slowly growing modes, 
while the theoretically fastest growing ones are too small for mean-field theory to be ap- 
plicable. As a result, we expect the dominant gravity- wave modes in the heat-salt system 
to be intermediate in size (20-30 fingers in width and 4-5 fingers in height), and grow 
with a non-dimensional growth rate of the order of about 0.04-0.05, taking into account 
that their real growth rate is expected to be smaller than the theoretically predicted one. 
Putting these results in dimensional form, using typical oceanic values oiv = 10~^m^s~"'^, 
KT = 1.4 X IQ'\t?s~^, ar = 2x IQ-'^K'^, g = lOms'^ and Tq^ = 10~2Km~\ we have 
d ^ 0.9 cm and IFGW — 8.2d, so that a typical finger is of the order of 7 cm wide and 
35 cm high. Hence the gravity waves expected to dominate the dynamics have horizontal 
and vertical wavelengths of the order of 1.5m - 2m, and grow on a timescale of the order 
of a few (3-5) hours. 

These gravity waves saturate rapidly as a result of localised breaking events, without 
directly inducing layer formation. Meanwhile, a spectrum of 7-modes also grow, with 
one of them destined to become the progenitor of the thermohaline staircase. Our results 
suggest that the dominant gravity-wave modes interact nonlinearly with the 7-modes 
and act as a filter for the smaller-scale ones. Only those with vertical wavelengths larger 
than that of the dominant gravity-wave modes are allowed to grow, but then appear 
to do so on a timescale well-predicted by mean-field theory. Since the fastest growing 
7-modes are those with the smallest vertical wavelength, we expect the progenitor of the 
staircase to be a mode with vertical spacing commensurate with the vertical wavelength 
of the dominant gravity-wave modes, i.e. around 5 fingers in height or about 1.5m-2m 
for parameters appropriate of the heat-salt system. This mode grows on a timescale of 
5-6 hours, and upon reaching a critical amplitude, causes regularly-spaced inversions in 
the mean density profile which are subject to direct overturning convection. A regular 
staircase then forms, with an initial layer spacing of about 2m. Note that once formed, 
experiments ( |KrishnamurtI||2003[ ) and theory ( |Radko||2003[ |2005[ |2007[ ) suggest that the 
staircase evolves further in time through successive mergers until an equilibrium layer 
depth is established. These events have not been observed in our 3D simulations yet since 
integrating the simulation on a merger timescale ( Radko|2005 ) is numerically prohibitive 
in 3D. 

Beyond giving new insight into the origin of oceanic staircases, our study marks a 
starting point in answering the question of whether gravity modes and 7-modes can 
be excited in other less readily observable fluid systems as well, and whether staircases 
may be expected or not. For example, it has been pointed out that vigorous fingering 
convection might occur in stars ( Vauclair|2004 Charbonnel fc Zahn|2007 Stancliffe et al. 



2007), where the resulting transport may account for observed peculiarities in luminosity 



or metallicity. Our work suggests that, despite the complexity of these highly nonlinear 
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processes, much can still be learned first by obtaining local flux laws for turbulent mixing 
in the parameter regimes typical of these systems, and then applying the mean- field 
theory framework as show in Paper I and in this work. 

A.T., P.G. and T.R. are supported by the National Science Foundation, NSF-093379, 
and T.R. is supported by an NSF CAREER. S.S. was supported by grants from the 
NASA Solar and Heliospheric Program (NNG05GG69G,NNG06GD44G,NNX07A2749). 
The simulations were partially run on the Pleiades supercomputer at UCSC, purchased 
using an NSF-MRI grant. Computing time was also provided by the John von Neuman 
Institute for Computing. We thank Gary Glatzmaier for many helpful discussions and 
for his continuous support. 

REFERENCES 

BoiS, P. A. & KuBlCKl, A. 2002 Double diffusive aspects of the convection in moist-saturated 
air. In Continuum Thermomechamcs , Solid Mectianics and its Applicattons , vol. 76, pp. 
29-42. Springer. 

Charbonnel, C. & Zahn, J. P. 2007 Thermolialine mixing: a physical mechanism governing 
the photospheric composition of low-mass giants. Astron. Astrophys. 467 (1). 

GuiLLOT, T. 1999 Interiors of Giant Planets Inside and Outside the Solar System. Science 
286 (5437), 72. 

Kraichnan, R. H. 1967 Inertial ranges in two-dimensional turbulence. Phys. Fluids 10, 1417- 
1423. 

Krishnamurti, R. 2003 Double-diffusive transport in laboratory thermohaline staircases. J. 

Fluid Mech. 483, 287-314. 
Krishnamurti, R. 2009 Heat, salt and momentum transport in a laboratory thermohaline 

staircase. Journal of Fluid Mechanics 638, 491 — h. 
Merceret, F.J. 1977 A possible manifestation of double diffusive convection in the atmosphere. 

Boundary Layer Meteorol. 11 (1), 121-123. 
Merryfield, W.J. 2000 Origin of ThermohaUne Staircases. J. Phys. Oceanogr. 30 (5), 1046- 

1068. 

OzGOKMEN, T., ESENKOV, O. & Olson, D. 1998 A numerical study of layer formation due to 
fingers in double-diffusive convection in a vertically-bounded domain. J. Mar. Res. 56 (2), 
463-487. 

Radko, T. 2003 A mechanism for layer formation in a double-diffusive fluid. J. Fluid Mech. 
497, 365-380. 

Radko, T. 2005 What determines the thickness of layers in a thermohaline staircase? J. Fluid 
Mech. 523, 79-98. 

Radko, T. 2007 Mechanics of merging events for a series of layers in a stratified turbulent fiuid. 

J. Fluid Mech. 577, 251-273. 
RUDDICK, B. & Gargett, A.E. 2003 Oceanic double-diffusion: Introduction. Prog. Oceanogr. 

56, 381-393. 

SCHMITT, R.W. 1979 The growth rate of super-critical salt fingers. Deep-Sea Res. 26A, 23-40. 
SCHMITT, R.W. 1994 Double Diffusion in Oceanography. Annu. Rev. Fluid Mech. 26 (1), 255- 
285. 

ScHMiTT, R.W., Ledwell, J.R., Montgomery, E.T., Polzin, K.L. & Toole, J.M. 2005 
Enhanced Diapycnal Mixing by Salt Fingers in the Thcrmocline of the Tropical Atlantic. 
Science 308 (5722), 685. 

Stancliffe, R.J., Glebbeek, E., Izzard, R.G. & Pols, O.R. 2007 Carbon-enhanced metal- 
poor stars and thermohaline mixing. Astron. Astrophys. 464, L57-L60. 

Stern, M.E. 1960 The salt fountain and thermohaline convection. Tellus 12 (2), 172-175. 

Stern, M.E. 1967 Lateral mixing of water masses. Deep-Sea Res. 14 (1), 747-753. 

Stern, M.E. 1969 Collective instability of salt fingers. J. Fluid Mech. 35. 

Stern, M.E., Radko, T. & Simeonov, J. 2001 Salt fingers in an unbounded thermocline. J. 
Mar. Res. 59 (3), 355-390. 



18 S. Stellmach, A. Traxler, P. Garaud, N. Brummell and T. Radko 

Stern, M.E. & Turner, J.S. 1969 Salt fingers and convecting layers. Deep-Sea Res. 16 (1), 
97-511. 

Tait, R.I. & Howe, M.R. 1971 Thcrmohalino staircase. Nature 231 (5299), 178-179. 
Tait, S. & Jaupart, C. 1989 Compositional convection in viscous melts. Nature 338 (6216), 
571 574. 

Toole, J. & Georgi, D. 1981 On the dynamics of double diffusively driven intrusions. Prog. 
Oceanogr. 10, 123-145. 

Vauclair, S. 2004 Metallic Fingers and Metallicity Excess in Exoplanets' Host Stars: The 

Accretion Hypothesis Revisited. Astrophys. J. 605 (2), 874-879. 
Veronis, G. 2007 Updated estimate of double diffusive fluxes in the C-SALT region. Deep-Sea 

Res. I 54. 

Walsh, D. & Ruddick, B. 1995 Double-diffusive interleaving: The influence of nonconstant 
diffusivities. J. Phys. Oceanogr. 25, 348-358. 



